Discrete spherical means of directional derivatives 
and Veronese maps 



Alexander Belyaev^, Boris Khesin'', Serge Tabachnikov 

Electrical, Electronic & Computer Engineering, Heriot-Watt University, Edinburgh, EHI4 4^S, UK 
''Department of Mathematics, University of Toronto, Toronto, ON M5S 2E4, Canada 
'^Department of Mathematics, Pennsylvania State University, University Park, PA 16801, USA 



Abstract 

We describe and study geometric properties of discrete circular and spherical means of 
directional derivatives of functions, as well as discrete approximations of higher order 
differential operators. For an arbitrary dimension we present a general construction for 
obtaining discrete spherical means of directional derivatives. The construction is based 
on using the Minkowski's existence theorem and Veronese maps. Approximating the di- 
rectional derivatives by appropriate finite differences allows one to obtain finite difference 
operators with good rotation invariance properties. In particular, we use discrete circu- 
lar and spherical means to derive discrete approximations of various linear and nonlinear 
first- and second-order differential operators, including discrete Laplacians. A practi- 
cal potential of our approach is demonstrated by considering applications to nonlinear 
filtering of digital images and surface curvature estimation. 



Introduction 

Mean value properties of functions play crucial role in analysis of many partial differ 



ential equations , numerical interpolation and integration 33 1 , computer tomography 
[23j . and many other areas of mathematics and engineering. In this paper we undertake 
a study of discrete circular and spherical means of directional derivatives of functions 
in view of their applications to the finite difference methods. We present a general con- 
struction for obtaining discrete spherical means of directional derivatives, analyze general 
properties of such means, and use them to derive finite difference approximations with 
good rotation-invariant properties for linear and nonlinear first- and second-order partial 
differential operators. The two principal components of our approach are based on the 
use of the Minkowski's existence theorem [sll and Veronese maps (l5l |. 

For the sake of simplicity, we focus on finite difference approximations of very basic 
differential operators, the gradient and Laplacian. We start from an elementary analysis 
of discrete circular means of derivatives and demonstrate how such discrete means can be 
used for practical estimation of surface curvatures (Section [T|) . Then, in any dimension, 
we obtain Minkowski-type formulas for general and regular grids by using the Veronese 
map (Section [2]). Finally, we analyze rotation invariance properties of discrete 2D gradi- 
ent, Laplacian, and related nonlinear operators defined on a square grid (Section[3]). 

Various finite difference approximations of the gradient and Laplacian are widely used 
in fiuid mechanics, electromagnetics, finance modeling, image processing, and many other 
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areas. Our interest in constructing reliable discrete approximations of these operators is 
threefold. Firstly, reliable and consistent approximations of the gradient and Laplacian 
are needed in various diffusion-type hydrodynamical models, such as the convection- 
diffusion equation 

ut + V ■ Vii = eAu , 

where the density u is convcctcd by a velocity field v and dissipates at different time 
scales. Its attractors exhibit very peculiar behavior as t — > oo and e ^ (see and 
we hope that the approximation scheme developed in this paper could complement the 
analytic treatment of the attractors. The Laplacian case is also particularly useful in 
computations related to the fast dynamo problem, where the proposed approximation 
formulas may simplify computations of the magnetic diffusion for iterations of a seed 
magnetic field, see 

Secondly, linear and nonlinear diffusion-type equations are widely used in modern 
image processing for a variety of tasks including edge detection, image restoration, and 
image decomposition into structure and and texture components j38ll2l|. As demonstrated 
in [39|, the use of discrete nonlinear diffusion operators with good rotation invariance 
properties is highly beneficial for such tasks. In addition, nonlinear image diffusion pro- 
vides us with a useful tcstbcd for the theoretical considerations below. Thirdly, accurate 
and reliable approximations of the surface gradient and Laplacian are key ingredients in 
many shape analysis studies. Below we discuss applications of the approximation formu- 
las both in image processing and geometry. Finally we note that the approach via the 
Veronese map described in the present paper can also be used for higher order differential 
operators, such as bi-Laplacian. 



1. Circular means of derivatives and their applications 

Continuous and discrete circular means. A simple way to demonstrate the rota- 
tional invariance of the 2D Laplacian 



consists of representing it as the circular mean (respectively, the spherical mean in 3D) 
of the second-order directional derivatives 
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The rotational invariance of the squared gradient |V/p can be demonstrated in a similar 
way as 
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More generally, the mean value representation for a directional derivative of / can be 
obtained from 
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where a and b are constants. Notice that ([3|) with a = d/dx and b = d/dy imphes ([T]) 
and setting a = df /dx and b = df /dy yields the equahty ([2]). 

Let us now derive discrete counterparts of the above mean value representations ([T}, 
([5]), and Consider a set of directions = (cos (^5^, sintpfc), fc = l,2,...,ri. Then a 
weighted sum of the second-order directional derivatives can be expressed as 



dxdy 
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Finding weights {wk} such that 

Y^WkC^'-^- = (5) 

yields, after a simple normalization {wk} — > {wk /J2'^k}, a discrete mean value rep- 
resentation for A/: 

and similar representations of the squared gradient 



and directional derivative 



+ b^ 
dx dy 



df 

Wh [a cos ipk + b sin ipk]^ , (8) 



respectively. Obviously the same set of weights {wk} can be used to obtain a discrete 
mean value representation for a quasi-Laplacian 

Note that the appearance of double angles in the formula exactly corresponds to the 
Veronese map V : (cos 0, sin (/)) i-^ (cos 2(^, sin 2(^), considered in Example 12.41 of the next 
section. 



Discrete mean value formulas and Minkowski problem for polygons. The key 

problem of finding normalized weights {wk \ J^'^k = 1} hi the relation ^WfeC^*'^'' = 
(see (O above) has a simple geometric meaning and is directly related to the Minkowski 
existence theorem (also called Minkowski's Problem, see e.g. [sH) for polyhedra. 
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Theorem 1.1 (Minkowski's existence theorem). Consider unit vectors ni,n2,...,nfe 

spanning R™ and a set of positive weights Wi, W2, . .. ,Wk- Then the equation 

wiUi + W2n2 + • • • + Wkiik = 

is a necessary and sufficient condition for existence of a convex polyhedron with facet 
unit normals rii , n2 , . . . , and corresponding facet areas Wi,W2, ■ ■ ■ ,Wk- Furthermore, 
this polyhedron is unique up to translation. 

The necessity part of this theorem is, of course, trivial and follows immediately from 
the Gauss divergence theorem. 

Let us note that in ([T]), and ^ one deals with the directions e{ip) = {cos (p, sin ip), 
< (/? < TT, which can be parameterized by the unit vector e'^^'^. Now consider a closed 
polygon such that the vectors e^*'^'= are the outward normals of the polygon's edges. Then, 
according to the Gauss divergence theorem, the edge lengths provide us with the desired 
set of weights {wk}- Vice versa, given a set of unit vectors {e^*'^'=} and positive weights 
{wk} satisfying the first condition in ([S]), Minkowski's existense theorem guarantees an 
existence of a convex polygon with edge lengths {wk} and outward normals e^*'^''. The 
left image of Figure [1] illustrates this geometric interpretation of ([5]) . 




Figure 1: Left: Minkowski's existence theorem for polygons provides us with a geometric interpretation 
of l(5j. Middle: a special case of a polygon circumscribed around the unit circle was considered in |21|1 in 
connection with the curvature estimation approach proposed there. Right: 3x3 stencils for estimating 
the gradient and Laplacian of functions defined on a square grid contain four directions and correspond 
to rectangles. 



In particular, if we consider a polygon circumscribed around the unit circle, as seen 
in the middle image of Figure (TJ the lengths of the polygon sides are given by 

tan/Jj_i +tan/3j, Pk = ipk+i - fk, 

and, therefore, the normalized weights are 

tan^j_i + tan^j 

Wj = "F^T -7, -TT-r- (10) 

^ ^(tan/3fc_i +tan/3fc) 

This set of weights was used in (2ll | for estimating the mean curvature on surfaces ap- 
proximated by dense triangle meshes. 

In Section [3] we show how to solve this key problem in any dimension and, in partic- 
ular, get a discretization of the corresponding Laplacian. 
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Application to curvature estimation. In the previous examples, we have considered 
function data defined on a regular grid. Here we deal with less organized data, namely 
with triangle meshes approximating smooth surfaces embedded into the Euclidean space 

Reliable estimation of curvature characteristics of polygonal surfaces is very impor- 
tant for many computer graphics and geometric modeling applications. Starting from a 
seminal paper of Taubin 1341 integral-based approaches are frequently used for curvature 
estimation purposes 21,2^221 (see also references therein) . Our approach to curvature 



estimation is based on using discrete circular means of directional surface derivatives. 
Given a smooth surface approximated by a dense triangle mesh, we employ such discrete 
circular means for estimating the mean curvature H = (fcmax + fcmin)/2 and the so-called 
curvedness R = \J (fc^iax + ^min) /■^i '^here fcmax and fcmin stand for the principal 
curvatures. (While the curvedness R is not a classical curvature measure, it is widely 
used in numerous applications including protein interaction analysis , heart physiology 
(40j . and human visual perception [35|.) 

Let n be a surface orientation normal. Denote by tmax and tmin the principal direc- 
tions corresponding to the principal curvatures fcmax and fcmin, respectively. For a point 
P on the surface, consider a unit tangent vector t((y5), which makes angle with tmax, 
and the corresponding directional curvature fc(<p). Similar to ([1} and ([5]) we have 
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respectively. Here (jlip immediately follows from Eulcr's formula 

k{ip) = fcmax COS (^^ -|- fcmin siu ip^ 

and p2p is obtained from a similar representation 



2 
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which, in turn, can be easily derived from the Rodrigues formulas 



dn dn 

~ fcmax tvi 



^, -^max ^max, ^. 

Ctmax ^^^min 

Similar to © and ([7]), let us consider discrete counterparts of pT|) and ([T^ 

H = J2^jHPj) and ^^ = Xl^j(a^^) ' (13) 

where weights {wj} satisfy ([5]) and are normalized by ~ 1. For the sake of sim- 

plicity, we use the set of weights (fTO)) proposed in [2l| . 

Consider now a triangle mesh approximating the surface and assume that each mesh 
vertex P is equipped with a unit normal n(P) which delivers a fairly good approximation 
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of the corresponding surface normal. Given a mesh vertex P and its 1-ring neighboring 
vertices {Qj}, as seen in the left image of Figure [2j simple approximations of tangential 
vectors at P and the corresponding directional curvature k{Lpj) and directional derivative 
of n arc given by 



PQj , 2PQ, ■ n dn_ n(Qj-) - n(P) 

IP -0,11' \\P-Q,\\'' dt, \\Q,-P\\ ' 



respectively. We refer to the right image of Figure [5] for a geometric explanation of 
the directional curvature approximation. Finally, (jl3p and (|14p provide us with discrete 
approximations of the mean curvatures H and curvedness R. 




Figure 2; Left: a mesh vertex P and its 1-ring neighboring vertices {Qj }. Right: a geometric explanation 
of the directional curvature approximation in I I14I1 . 



FigurcOprcscnts color encodings for the mean curvature H and curvedness R detected 
on a geometrical model with many small-scale details (see the electronic version of this 
paper for a color version of this figure) . 




Figure 3: Color encodings (please see the online version of this paper for this figure in color) for the 
mean curvature (left) and curvedness (right) detected on Gargoyle model (center). 
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A visual comparison of (fT3|) and (jl4[) with the curvature estimation methods devel- 
oped in [s^l and [2^ indicates that our discretization is less robust to noise but better 
reflects small-scale surface details. One can improve the approximation by using more 
sophisticated approximations for the tangent directions, directional curvatures, and di- 
rectional derivatives of the normal instead of common-sense estimates (|14p . Here our goal 
was to reveal potential applications of our approach to curvature-based shape analysis. 

2. Discrete spherical means of derivatives 

General construction and Veronese maps. In this section we describe a general- 
ization of the discrete mean value formulas of the preceding section from the circle case 
to arbitrary higher dimensions. In particular, we prove the following 

Theorem 2.1. Given a degree fc 7^ and a collection of points Pi,i — l,...,N, in 
general position on the unit sphere in C R""*'"'^, with N > {"'^'') — ("^^^^); there exist 
weights wi (not all equal to zero), which can be found constructively, such that for every 
harmonic homogeneous polynomial F in R""*"^ of degree k one has '^^Wi ■ F{Pi) ~ 0. 
If N > {"'^'') , the same statement holds for any homogeneous polynomial F in R"''"^ of 
degree k. 

For fc = 1, we have N > n, and this theorem can be viewed as a direct corollary 
of the Minkowski theorem. Note that for a possible application to computations of the 
mean curvature of a manifold, the function F can be the sum of a constant {k = 0) and 
a homogeneous quadratic term (fc = 2), and one needs to find the weights for obtaining 
the mean value F of such a function from its values at the given set of points. (Of course, 
all nonzero harmonics should give us the zero mean.) The cases of computations for the 
discrete Laplacian and square gradient operators are similar. 

For instance, we use Theorem l2.1l to construct a discrete Laplacian as follows. 

Definition 2.2. Let O be a vertex in an arbitrary irregular grid with neighbouring 
vertices Pi . Define the discrete Laplacian of a function / at a point O to be 

iA/(0):=^u;,|^, (15) 

where points Pi are normalized from Pi to lie on the unit sphere centered at O, the weights 
Wi are given by the above theorem for fc = 2, and d'^f/de^, are discrete approximations 
of the second derivatives of / at O in the direction OPi . 

Definition (|15p is inspired by formula ^ and, as we will see later, generalizes the 
standard discrete Laplacians defined on regular grids. 

Remark 2.3. The result of Theorem 12.11 is very close to constructions of so called 
"spherical designs." The latter is a finite set of points on S*" such that the average value 
of any polynomial F of degree less than or equal to k on this set equals the average value 
of the polynomial F on the sphere. The concept of a spherical design was introduced in 
lol . The existence and structure of spherical designs in a circle was studied in [l^ ■ In 
32| the existence of designs of all sufficiently large sizes was proved: there is a number 
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N{n,k) depending on n and k, such that a design exists for any N > N{n,k). Many 
examples and recent results can be found in the survey Q- 

The statement of Theorem 12.11 is somewhat different: we do not require the weights 
to be equal. In many applications, however, one uses the regular latices, and the weights 
turn out to be equal, as we discuss below. 

Proof— Construction. First, we prove Theorem 12.11 for a linear functional L (i.e. for 
k = 1). In this case the solution is given by the Minkowski theorem ll.il Namely, consider 
the tangent hyperplanes to the sphere S"" at points Pi; they form a polyhedron. Here 
we use the fact that Pi,i = 1, . . . , N, are in general position, so that no hyperplanes are 
parallel and no more than n + 1 intersect at a point (note that the polyhedron is not 
necessarily circumscribed about the sphere, but the hyperplanes containing its facets are 
tangent to the sphere). Let Wi be the signed codimension-one volumes of the faces of 
this polyhedron. Then Wi ■ Pi = according to the Minkowski theorem (where we 
understand Pi as vectors in R""*"^ and where the weights may be negative), and hence 
J2i Wi • L{Pi) = for every linear function. 

Next, consider harmonic polynomials of degree k in R"+^. Their restrictions on the 
sphere form the space Jik of spherical harmonics of degree k. Consider the harmonic 
Veronese map V : R"+^ R', where V{xi, Xn+i) = {xi — x|,...), that is, the 
map defined by taking a basis of harmonic polynomials in xi, ...,Xn+i of degree k. The 
dimension is given by the formula: q = dimHk = ("^'') - Then F = L oV 

where L is a linear function on R"^. Let Wj be the constructed above weights for the 
collection of points V{Pj) in R"^. This gives us the required formula Wi ■ F{Pi) ~ 0. 
Note that the general position condition on Pi now means that the points Qi = V{Pi) 
in R'' , after normalization, satisfy the above condition that the tangent planes to the 
sphere 5''^^ form a polyhedron. 

Finally consider all homogeneous polynomials in R"+^ of degree k. Their restrictions 

on the sphere form the space T-Lk®T-Lk-2®T-Lk-i® The dimension of this space is {^^^) , 

and we repeat the same construction involving the Veronese map for each homogeneous 
component as before. 

Example 2.4. Let us demonstrate this approach in the case of a circle and quadratic 
functions. The harmonic Veronese map for quadratic polynomials sends the coordi- 
nate functions {xi, . . . , x„+i} to the space of quadratic harmonic polynomials (whose 
restrictions on the sphere give us spherical harmonics). For the circle C R^ 
this harmonic Veronese map is V : R^ — )■ R^, where {x,y) i— > (x^ — y^,2xy), or 
(cos sin (/)) i-> (cos 20, sin 2(/)), see a detailed analysis in the preceding section. 

Note that if we map R^ ^ to the full 3-dimensional space R^ of quadratic polynomials 
{x"^ , , 2xy) , then the image of the circle C R^ will lie in the plane R^ C R'^(due to 
the relation x^ + y"^ = 1), so the point images V{Pi) will not be in general position to 
yield a polyhedron circumscribed around a 2-sphere in R^. 

Example 2.5. For the case of 2-sphere and quadratic functions the space of harmonic 
polynomials is 5-dimensional, and our harmonic Veronese map sends S"^ to this R^. The 
point images lie in general position, and sufficiently many points (at least 6) would yield 
a polyhedron, which allows one to find the corresponding weights Wi. 

Remark 2.6. Since the points Qj — V{Pj) should lie on the sphere S''~^ , rather than be 
generic vectors in R'^, one can project them to the sphere by rescaling. Namely, let Qj be 
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a (sufficiently generic) collection of non-zero vectors in R"?. Assign the following weights 
Wj to these vectors. Consider the collection of points Qj/\Qj\ on the unit sphere 5"^"^, 
draw the tangent hyperplanes to the sphere, and let uj be the signed (g— l)-dimcnsional 
volumes of the faces of the resulting polyhedron. 

Proposition 2.7. For every linear function L on , one has: Wj-L(Qj) = 0, where 
the weights are Wj := Uj/\Qj\. 

Proof. The Minkowski formula says that X^^j ' Qj/\Qj\ — 0, hence = J^"^] ' 
LiQ,)/\Q,\=j:^Wj-L{Q,). QED 



Weights for regular lattices. Now, from the general consideration of points in general 
position we move to regular ones, and consider the set of all "m- mid-points," i.e., the set 
Mm C S" of points whose coordinates are all possible combinations of m nonzero coordi- 
nates equal to ±1 among n spots, starting with (±1, ±1, 0, 0) to (0, 0, ±1, ±1), 
which arc scaled to belong to the unit sphere S*" C R""^^. Here m is any fixed integer 
between 1 and n. We also confine to the case of quadratic harmonics, fc = 2 in the 
general problem above. 

Example 2.8. In dimension 3, for m = 1 we get Mi consisting of 6 midpoints of 
the cube faces, which are also the vertices of an octahedron: (±1, 0, 0), (0, ±1, 0) and 
(0, 0, ±1). The set M2 consists of 12 scaled mid-points of the cube edges: (±1, ±1, 0)/^/2, 
(±1, 0, ±l)/-\/2, and (0, ±1, ±l)/y/2. The set M3 contains exactly 8 vertices of the cube: 
(±l,±l,±l)/^/3. 

Theorem 2.9. Let F{xi, . . . ,Xn+i) be a homogeneous quadratic function, M„i C 5" be 
the set of m-mid-points for any fixed m, I < m < n. Then 

where vol{S^^) is the volume of n- dimensional sphere and ^Mm is the number of points 



One can see that in all these "regular cases" one counts the points of Mm with equal 



weights. This theorem generalizes the corresponding 3D consideration (see e.g. j24l|). 
which gives three basis stencils used in the Laplacian computations. 
Proof. Write F = c + G, where c e R and G a harmonic quadratic polynomial. For a 
constant c the claim is obvious. For G one has Jg„ GdS — 0, so we need to show that 
12xeM ^i^) ~ ^ every harmonic quadratic polynomial G. 
Consider the full Veronese map V : R" — > R'^+^, given by 

(Xi, ...,Xn+l) {xl, ...,x'^^i;XiX2, ...,XnXn+l) , 

with q + 1 = (77- + l)(n + 2)/2. Let H be the space of linear functions on R''+^ that vanish 
on the vector ^ — (1,...,1;0,...,0) with first ri + l nonzero components (we separate them 
by semicolon). Then, for each quadratic harmonic G, one has G = HoV for some H € T-l. 
Indeed, harmonic polynomials form a hyperplane in the space of linear functions on R'^+^ 
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being given by one linear relation (coming from the sphere equation + ... + x'^_^_l = 1), 
and the polynomials — X2, ■■■,x'fi — x'^^j^,xiX2, ...jXnXn+i are harmonic and equal to 
on this vector ^. Thus we need to check that X]QGy(M ) H{Q) — ioT all H E H and 
any m = I, ...,n + 1. 

For instance, for the set Mi, which consists of (0, 0, ±1, 0, 0)-typc points, the 
image V{Mi) consists of vectors Qi: (0, 0, 1, 0, ...0; 0, ..0) with I's in one of the first n 
coordinates. In the case of AI2 with 2n{n — 1) "diagonal" points 

(0,...,0,±l,0,...,0,±l,0,...,0)/^/2, 

the image V{M2) consists of vectors Qi: 

(0, ...,0, 1,0, ...0, 1,0, ...,0 ;0, ...,0,±1,0, ...,0)/2 

with all possible pairs of I's among first n + 1 coordinates and the ±1 at the respective, 
ij-place. Similarly, for any m one can see that in each of these cases we have ^^^'^^ = 
const ■ C Hence H{Qi) = for aU H e H, as claimed. QED. 

Thus, while Theorem 12.11 allows one to find the weights Wi for a generic set of vec- 
tors Pi, Theorem 12.91 gives the explicit weights in the case when F is a homogeneous 
quadratic polynomial (i.e., fc = 2) and vectors Pi correspond to the m-mid-points of the 
n-dimensional cube, 1 < m < n. 

Remark 2.10. For computations we are often interested mostly in the case C R^. 
Here is the explicit statement for M = M2 adjusted to this case: let F{x, y, z) be 
a homogeneous quadratic function, M C S"^ be the set of 12 points: (±1, ±1, 0)/a/2, 
(±1,0, ±1)772, and (0,±1,±1)/V2. Then 

-L / F dS^^y F{P). 
47r 12 ^ ^ ' 

Thus in all of the above cases we need to use the same weights on these points. 
This equality of weights for a given set Mm can also be obtained from the symmetry 
consideration. One can also take combination of the above lattices Mm with with different 
weights for different m (provided that the total sum of the weights is equal to 1), e.g., 
by considering arbitrary combinations of the 3 lattices in R"^: mid-edges, vertices of the 
cube, and at the middle of faces. By imposing additional requirements, one can optimize 
these coefficients: for instance, to achieve an isotropic approximation, or with the next 
approximation smallest in an appropriate sense. 

Remark 2.11. Points of V{Mi), evidently, form the coordinate [n + l)-simplex in the 
space spanned by first n± 1 coordinates, which is a slim subspace in R^+^. In the case of 
M2, the points of V{M2) form the standard q-simplex in a hypcrplane in R^+^. Indeed, 
in the latter case, we see that for all Qi G V{M2) one has [Qi, Qj) = 1/2 for i 7^ j and 
{Qi,Qi) = 3/4, i.e., the vectors Qi have the same length and same angles between them. 

Remark 2.12. Above we considered two discretization procedures, both using the 
Veronese map: a generic set of points {Pi} and the mid-points sets. While for a generic 
set in R'?+^ the coefficients are given by the Minkowski formula for the corresponding 
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polyhedron in the image, it is not the case in the mid-point cases. Indeed, due to a very 
regular structure of the image points Qj = V{Pj), the corresponding tangent planes go 
to infinity, and the corresponding volumes of the faces are infinite. 

In a sense, if the image points Qj are linearly dependent in R''+^( mod ^), and the 
corresponding polyhedron has faces going to infinity (like in the standard basis case) , we 
have to quotient out these "directions to infinity," and consider the volume of faces in a 
polyhedron of smaller dimension. 

Example 2.13. For n = 2 and for the standard basis we did not get a polygon on 
the plane = R^( mod ^). Instead, it reduces to a pair of points in R^. It is not 
enough to start with two generic points in this case, but if we take three generic points 
on the circle, we obtain an interpolation formula. This is a manifestation of the general 
phenomenon described in the following proposition. 

Proposition 2.14. There is no discrete mean value formula for quadratic harmonics 
and a given generic collection of n + 1 points on S" . 

Proof. Let Pi,..., P„+i £ 5", and let V : R"+i ^ R^+i with q + l = (n + l){n + 2)/2 
be the full Veronese map: 

Let F(Pfc) = Qk. Let ^ = (1, 1; 0, 0) S R«+i. 

We want to find constants Wk such that ^ ■ Qk = This is equivalent to the 
following equality: 

Dtagiwk) = iA*A)-\ 

where A is the matrix made of the vectors Pi, P„_|.i, as the following linear algebra 
shows. Indeed, the equality J^'^kQk = ^ equivalent to 

^WiXk,iXi.i = Sk,i, 

i 

which is equivalent to A ■ Diag{w) ■ A* — E, the identity matrix. This is equivalent to 
{A*A)~^ = Diag{w), as claimed. 

Of course, if {Pk} is an orthonormal frame then A is orthogonal, and all iffc = 1, as in 
the theorem above. The matrix AA* is diagonal if and only if the vectors Pi are pairwise 
orthogonal. This is the case for the set Mi discussed above, which consists of the basis 
points, but, of course, not true for generic collection of 71 + 1 points on S". This implies 
that one cannot have an interpolation formula with n + 1 generic points. QED 

Note that on the other hand, if one has sufficiently many points, so that their Veronese 
images are linearly dependent in j{,"("+i)/2 ^ mod ^), then any such linear relation pro- 
vides the weights that yield an interpolation formula. 

Higher-dimensional Laplacian approximation formulas. The results of the present 
section allow one to easily construct various discrete Laplacians and squared gradient 
similar to representations ^ and ([T]), respectively. In the beginning of the section we 
discussed the construction of the discrete Laplacian for an arbitrary grid, see formula 
03: 
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To obtain the weights Wi for a point O with neigboring vertices Pi we first normalize the 
latter to get Pi lying on the unit sphere centered at O. Then we consider the Veronese 
map and obtain points Qi. Finally we obtain the weights Wi — Ui/\Qi\ by finding the 
signed (g — l)-dimensional volumes Ui of the faces of the polyhedron in the Q-space, 
see Proposition 12.71 An analog of the formula ([7]) for the square gradient instead of the 
Laplacian is completely analogous. 

For a regular grid, Theorem 12.91 delivers similar formulas for a discrete Laplace oper- 
ator where all weights Wi are equal: 

A/(0) = C™^|^, (16) 

where points Pi belong to the grid of m-midpoints, Pi S M^, and the constant Cm 
normalizes the sum, Cm = #-^i/#-Wm, by comparing the number of neighbours in Mm 
with that in the coordinate grid Mi, the most straightforward definition of the discrete 
Laplacian. 

Furthermore, one can consider linear combination of the Laplace approximations (jl6p 
with different coefficients a™ satisfying X)m=i '^m ~ 1- -^^Y such linear combination 
delivers a different discrete Laplacian and one can run a separate optimization problem on 
these coefficients am to find a discrete Laplacian with better rotation-invariant properties. 
We discuss the work in this direction for the 2D Laplacian in the next section. 



3. Discrete circular means of derivatives and isotropic finite differences 

In this section we show how to use discrete circular means of Section [1] for approx- 
imation of various differential operators in 2D by finite differences with good rotation 
invariancc properties. 

The simplest way to construct such approximations consists of substituting in ([5]), ([7]), 
([5]), and (O three-point finite differences instead of the first- and second-order derivatives. 
For example, constructing a discrete Laplacian corresponding to (jG]) can be done as 
follows. Consider a point p inside a convex bounded domain il. Let the straight line 
determined by p and direction intersect the boundary dVl in two points qi and q2. 
Denote by pi and p2 the distances from p to qi and q2, respectively. The second-order 
ej^-directional derivative of /(p) can be approximated by 



del 



PlP2 
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(17) 



where pi and p2 are assumed to be small. Now a discrete approximation of the Laplacian 
is obtained by summing up (|17l) over the set of directions e^^, k — 1, 2, . . . , n. 

Discrete Laplacian ([5]), PT)) can be considered as a special case of a general construc- 
tion studied in [s^ (see also references therein). 

It is interesting to note that integration of (|T7)) with respect to (p 
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(18) 



leads to an extension of a pseudo-harmonic interpolation scheme proposed in [13|. It 
can be shown that, given /(q) defined on dil, (fTS)) and its generalizations deliver the 
harmonic interpolation if fi is a circle (a ball in the multidimensional case) 

Unfortunately combining circular means with directional finite differences docs not 
automatically guarantee rotation invariance properties of the approximating differential 
operators. The reason is that the finite differences introduce directional errors. However, 
as seen below, for a regular grid, such directional errors can be uniformly distributed 
over the directions by choosing appropriate weights in discrete circular means ((5]), (O, 
(0), and dH). 

Gradient and Laplacian stencils for square grids. Now let us focus on analyzing 
rotation invariance properties of regular grid stencils corresponding to ([5]) and ([7]). 

Consider a two-dimensional square grid with step-size h ^ 1. Each grid vertex has 
eight nearest neighbors: two horizontal, two vertical, and four diagonal. Thus, according 
to our geometric interpretation of ([5]), ©, and ([5]), we are given four directions and 
corresponding unit normals 

{e^'"^}, 9? = 0,7r/4,7r/2,37r/4 

which define a rectangle aligned with the coordinate axes. 

Consider such a rectangle with the edge lengths 2 and w, as shown in the right image 
of Figure [1] and replace the first- and second-order directional derivatives in ([B]) , ([7]) , 
and (HI by corresponding central differences. We arrive at the following stencils for the 
x-derivative and Laplacian 



d_ 

dx 



2h{w + 2) 
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(19) 



and 



1 w 1 

w — 4(w 4-1) w 



1 



w 



1 



(20) 



respectively. The stencil Dy for d/dy is obtained from Dx by 7r/2-rotation. Here and 
everywhere below the 3x3 matrices are understood as finite difference operators acting 
on functions defined on a square grid with step-size h. 

Formulas ([Tn)) and (PH)) provide us with consistent parameterizations of the 3x3 
stencils for the gradient and Laplacian. 

Remark 3.1. The standard central difference and five-point Laplacian is obtained 
for w = oo. Setting w = 1 leads to the so-called Prewitt masks for the Ist-order 
derivatives and Laplacian 2^, 12 1 . The case w = 2 is commonly used in image processing 
applications and yields the standard Sobel mask for the derivative [2^, [l2| and a 9-point 
discrete Laplacian, which possesses good isotropic properties [l3|. The case w = 4 was 
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analyzed in where it was shown that 



w=4 



1 

~8h 



-1 1 
-4 4 
-1 1 



d 
dx 



4 
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(21) 



(22) 



as /i — >■ 0. Since the Laplacian is isotropic, the right hand sides of (PT|) and ([^^ de- 
hver asymptotically optimal asymptotically rotation-equivariant 3x3 stencils for the 
x-derivative and Laplacian, respectively, as <C 1. In particular, this explains why ((20)) 
with w = 4 is often used for a numerical solution of the Poisson equation (loL Chapter 
3, §1], @, Chapter 3, §10]. 

Discrete nine-point Laplacian (j20p can be represented as a linear combination of two 
basis five-point Laplacians 



= oiL^ + /3Lx with a 



w + 2 



1 

7^ 



1 

-4 1 
1 



and 



and f3 = 
1 

^ 2^ 



where 



Setting w — 4 yields a = 2/3 and /3 = 1/3. 

In Figure m we use a Gaussian bump function 

g{x,y) = exp (x^ + y^) 

to test isotropic properties of four discrete stencils 



(23) 



, and 

W — 2 10—4: 

. As (|22|) suggests, with w = 4 



Namely, the figure displays contour lines of (A — L^) 
demonstrates the best performance if the grid step-size h is sufficiently small. 

Similar results arc obtained for the x-derivative 3x3 stencils (jl9p and corresponding 
y-derivativc stencil. Analogously to the case of the Laplacian, expansion (j2ip implies 
that p9|l with w = 4 is rotationally optimal for sufficiently small h. 



Remark 3.2. It is easy to increase the approximation accuracy of (PT|) and 
example, one can rewrite (PT|) as 



For 



— - (l —A 

dx V 



0{h^) 



1^2 ; ^\w=4 



and use stencil (|22|) for approximating A. 
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Quasi-Laplacian stencil with good rotation invariance properties. Consider now 
a quasi-Laplacian operator 

L" = V • [a{x, y)V] = aA + Va • V (24) 

where a{x,y) is given. In view of ([22]), dH]), and the optimal 3x3 stencils for 
are obtained from and when it; = 4. 

In practice, the diffusion coefficient a{x,y) in is defined at the staggered grid 
points which lie on halfway between the standard grid points. Following [H, A. 3. 2] let us 
approximate (|24p using a linear combination of two basis stencils 

all + ^Ll 

with 

= 7^ ["+o/,+i,, +a_o/.-i,, +ao+/,,3+i +ao-/.,,-i 

and 

Lxif] = ^ +a-+/.-i,,+i +a+-/.+i.3-i +a-./>-i,.-i 

where a^„, a^^, and a^^ denote the values of a(x, ?/) at {{i ±l/2)h,jh), {ih,{j ±l/2)h), 
and ((i ± l/2)h, {j ± l/2)/i), respectively. 

Assume h <^ 1. Since the same weights are used in the discrete directional mean value 
representations ([6]) and ([9]) we can expect that the optimal rotation-invariant stencil for 
quasi-Laplacian is obtained when a = 2/3 and /3 = 1/3, or equivalently, when w = 4. 
This indeed turns out to be the case in the lower order terms: 

Proposition 3.3. The discrete approximation aL'^ + f3L'^ of the quasi-Laplacian L" is 
is asymptotically rotation- equivariant for a = 2/3 and 13 = 1/3 up to order O (/i'^) as 
h->-0. 
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Proof. A straightforward calculation shows that 



L-[f] = V ■[aix,y)yf] 
+ T2 



D,, ^ (a, /) + 2D, , (a, /) + (a, /) + ^D, ^ (a, /) 



(25) 



with 

-02,2(0,/) 



9a 9^/ ^ 9a 9^/ 



dx^dy^ dy^ 

da d^f da d^f 



dx dx^ dy dx'^dy dx dxdy^ dy dy^ ' 
d^a 92/ d^a 92/ 4 d^a 9^/ 1 92a 9^/ 1 d'^a 9^/ 



9x2 ^2:2 Qy^ Qy2 

d^a 9/ 9^a 9/ 



3 dxdy dxdy 3 9a;2 9j/2 
9^a 9/ 9^a 9/ 



3 9y2 9^2 



9a;'^ dx dx^dy dy dxdy^ dx dy^ dy 



It is interesting that expansion (|25l) does not contain ft.'^-terms. 

Let us now demonstrate that these differential quantities are rotationally invariant. 
Obviously Dg^{a,f) equals a{x,y)A^f which is rotationally invariant. Let us consider 



da d^f 



de^ del 



d^a 9V 
9e2 del 



^3,1 (a,/, 



d^a df 



del '^^v 



Note that -^{a, f,(p) = R,,^{f,a,ip) and D,,{a,f) — D,,{f,a). Direct calculations 
show that 



-Ri,3(«:/' dip 



T 



D,,{aJ) and / 
Jo 



,{a,f, (p) dip 



^2,2 («:/)■ 



This completes our proof that the h^-tcim in ([25]) is rotation invariant. QED 

Application to nonlinear diffusion image filtering. Accurate numerical implemen- 
tations of nonlinear diffusion processes governed by 



du{x,y,t) 
dt 



div {a{x, y, u, Vu)Vu) , u{x, y, 0) = uo{x, y), 



(26) 



subject to appropriate boundary conditions, are important for a number of disciplines 
including computational physics, magnetohydrodynamics, financial mathematics, and 
image processing. In particular, in image processing, adaptive image smoothing is often 
carried out by ([^^ with 



a(x,?/,u, Vu) = exp(— I Vu|/A), 



(27) 



as suggested by Perona and Malik in their seminal paper [23. Figure [5] demonstrates how 
our numerical implementation of (j26p and (j27p performs an adaptive image smoothing 
and removes texture from an image. 
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Figure 5: Demonstrating the power of nonlinear diffusion for image filtering. Left: original image "trui". 
Right: after filtering with II26II and II27II . Middle: a comparison of onc-dimcnsional image slices of the 
original and filtered images. 



In image processing, the system ([26 | . ((27)) is typically solved by finite differences meth- 
ods with explicit or implicit Euler schemes and several dozens of iterations are usually 
required in order to achieve desired image filtering effects. For a number of applications 
including image segmentation and feature extraction, accurate estimation of the gradi- 
ent directions in an image filtered by (j26p and (|27[) is required. Unfortunately standard 
finite difference schemes have poor rotation-invariance properties. If such schemes are 
used repeatedly, the directional error is accumulated. See 37, 3^ for examples of finite 
difference schemes with improved rotation invariance properties. 




Figure 6: Left: the contour lines of Gaussian bump function II23II arc perfect circles. Middle: contours 
of the Gaussian filtered using our implementation of nonlinear diffusion I I26I I. I I27I I. Right: contours of 
the Gaussian filtered using code from [2^ . 



Our approach to numerical solving of (j26p . (|27p combines the explicit Euler scheme 
with obvious generalizations of the optimal rotation-invariant 3x3 stencils derived in pre- 
vious sections for gradient, Laplacian, and quasi-Laplacian. To demonstrate advantages 
of our approach we apply our numerical implementation of (^5)) . (P?]) to (^5]) . Figure [S] 
shows level sets of before and after applying nonlinear diffusion (^5]) . (P?]) . We also 
compare our implementation with that from from (26| , where five-point stencils are used 
for numerical solving (^5]) and (P7|) . 

Frequency response analysis. Now we take a brief look at stencils (|19l) and (PH)) from 
a frequency point of view. This standpoint is important when dealing with phenomena 
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possessing a range of space and time scales. Turbulent fluid flows, electromagnetic and 
acoustic wave propagation, multiresolution analysis of financial data, and multiscale im- 
age processing are among common examples. 

For the sake of simplicity assume that our square grid has a unit step-size: h = 1 
(this can be easily achieved by rescaling) . One can see that (fTO|) consists of the simplest 
central difference operator combined with a smoothing kernel applied to the orthogonal 
direction: 

Dx = ^\ -I ll-^--7[l w iV. (28) 

The eigenvalue (also called the frequency response) of ([T^ corresponding to the eigen- 
function exp(i(wia; -I- U!2y) is given by 

w + 2 



I smwi 



2 cos U!2 



-TT < Wfe < TT, k = 1, 2, 



where the first term of the product is the frequency response for the central difference 
and the second one corresponds to the smoothing kernel in ([25)) . 

It is well known that the frequency response i sin w delivers a satisfactory approxi- 
mation of iuj, the frequency response for the ideal derivative, only for sufficiently small 
frequencies w (see, for example, (13 . Section 6.4]). Now it is clear how ([T9| improves the 
central difference: smoothing due to the use the central difference operator instead of 
the true x-derivative is compensated by adding a certain amount of smoothing in the 
y-dircction. Thus (jl9|) and its y-direction counterpart do a better job in estimating the 
gradient direction than in estimating the gradient magnitude. 

If the goal is to achieve an accurate estimation of both the gradient direction and 
magnitude, we can combine (jl9|) and (|20)) as follows. Let 



S = 



be the 3x3 identity kernel. Note that 
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which can be considered as simultaneous smoothing (averaging) with respect to both the 
coordinate directions. Thus it is natural to use 



1 



w + 2 



(29) 



which combines (jl9|) with a Laplacian-based sharpening. The frequency response func- 
tion corresponding to ([2^ is given by 



H(uj) ~ i sinoj 



w + 2 



w + 2 cos uj 



(30) 



where parameter w should be chosen in such a way that (PH)) delivers a good approxi- 
mation of i UJ, the frequency response for the true derivative. For example, the Taylor 
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series expansion of ([30)) with respect to a; at a; = shows that the best approximation 
of the exact x-derivative for w ^ 1 is achieved when w = 4 which corresponds to a Pade 
approximation. However this choice of parameter w is not optimal if wc deal with a 
wider range of frequencies. 

Note that (|29p is equivalent to an implicit finite difference scheme 



1 



w - 



2 



(31) 



Implicit finite differences are widely used for accurate numerical simulations of physical 
problems involving linear and non-linear wave propagation phenomena, see, for example, 
0, [23, Section 5.8]. An implicit finite difference scheme is usually characterized by 
its frequency- resolving efficiency, the range of frequencies uj over which a satisfactory 
approximation of the corresponding differential operator is achieved. 

Motivated by 3^ 17 1 where a frequency-based analysis of simple 3x3 stencils for 
first-order derivatives was performed, we found out that ([501 ([HT|) with w = 10/3 deliver 
a good frequency resolving efficiency for a sufficiently large frequency range, as seen in 
Figure [71 Of course, the problem of finding an optimal w is application-dependent and 
the optimal value of w depends on optimization criteria used and a range of frequencies 
where the optimization is sought. 



3_ 



■ exact differentiation 

- w = CO (2nd-order central difference) 

■ w = 4 
■w = 10/3 




Figure 7: Visual frequency-response analysis of II29II for various values of parameter w. 
frequency-resolving efficiency of II30I I with w = 10/3. 



Notice good 



The approach via the Minkowski theorem and Veronese map developed in Section [2l 
is suggestive for an extension of the above results to 3D and to grids in any dimensions. 
We leave this for future publication. 

Acknov^fledgments. B.K. is grateful to the Ecole Polytechnique in Paris and the 
Max Planck Institute in Bonn for their hospitality during completion of this paper. He 
was partially supported by an NSERC research grant. S. T. was partially supported by 
the Simons Foundation grant No 209361 and by the NSF grant DMS-1105442. 

19 



References 



[1] Arnold, V., Khesin, B., 1998. Topological methods in hydrodynamics. Springer. 
[2] Aubert, G., Kornprobst, P., 2006. Mathematical Problems in Image Processing, 2nd Edition. 
Springer. 

[3] Bannai, E., Bannai, E., 2009. A survey on spherical designs and algebraic combinatorics on spheres. 
European Journal of Combinatorics 30 (6), 1392-1425. 

[4] Belyaev, A., June 2006. On transfinite barycentric coordinates. In: Proceedings of Fourth Euro- 
graphics Symposium on Geometry Processing. Sardinia, Italy, pp. 89—99. 

[5] Bickley, W. G., 1948. Finite difference formulae for the square lattice. Quart. J. Mech. Appl. Math. 
1, 35-42. 

[6] Birkhoff, G., Lynch, R. E., 1984. Numerical Solution of Elliptic Problems. SIAM. 

[7] Bradford, J. R., Westhead, D. R., 2005. Improved prediction of protein-protein binding sites using 
a support vector machines approach. Bioinformatics 21 (8), 1487-1494. 

[8] Childress, S., 1992. Fast dynamo theory. In: Topological aspects of the dynamics of fluids. Kluwer 
Academic Publishers, Dordrecht, pp. 111-147. 

[9] Colonius, T., Lele, S. K., August 2004. Computational aeroacoustics: progress on nonlinear prob- 
lems of sound generation. Progress in Aerospace Sciences 40 (6), 345—416. 
[10] Delsarte, P., Goethals, J. M., Seidel, J. J., 1977. Spherical codes and designs. Geometriae Dedicata 
6, 363-388. 

[11] Friedman, A., Liftman, W., 1962. Functions satisfying the mean value property. Transactions of 

the American Mathematical Society 102 (1), 167-180. 
[12] Gonzalez, R. C, Woods, R. E., 2008. Digital Image Processing, 3rd Edition. Pearson Prentice Hall. 
[13] Gordon, W., Wixom, J., 1974. Pseudo-harmonic interpolation on convex domains. SIAM J. Numer. 

Anal. 11 (5), 909-933. 
[14] Hamming, R. W., 1998. Digital Filters, 3rd Edition. Dover. 
[15] Harris, J., 1992. Algebraic Geometry: A First Course. Springer. 

[16] Hong, Y., 1982. On spherical t-designs in r^. European Journal of Combinatorics 3, 255—258. 

[17] Jahne, B., Scharr, H., Korkel, S., 1999. Principles of filter design. In: Handbook of Computer Vision 

and Applications, Vol. 2, Signal Processing and Applications. Academic Press, pp. 125—151. 
[18] Kamgar-Parsi, B., Kamgar-Parsi, B., Rosenfeld, A., 1999. Optimally isotropic Laplacian operator. 

IEEE Transactions on Image Processing 8 (10), 1467-1472. 
[19] Kantorovich, L. V., Krylov, V. I., 1956. Approximate Methods of Higher Analysis. Noordhoff- 

Interscience. 

[20] Koenderink, J. J., Van Doorn, A. J., 1992. Surface shape and curvature scales. Image and Vision 
Computing 10, 557-565. 

[21] Langer, T., Belyaev, A., Seidel, H.-P., 2007. Exact and interpolatory quadratures for curvature 
tensor estimation. Computer Aided Geometric Design 24 (8-9), 443—463. 

[22] Merigot, Q., Ovsjanikov, M., Guibas, L., 2009. Robust Voronoi-based curvature and feature es- 
timation. In: Proc. of SIAM/ACM Joint Conference on Geometric and Physical Modeling, pp. 
1-12. 

[23] Natterer, F., 2001. The Mathematics of Computerized Tomography. SIAM. 

[24] Patra, M., Karttunen, M., 2005. Stencils with isotropic discretisation error for differential operators. 
Numerical Methods for Partial Differential Equations 22 (4), 936-953. 

[25] Perona, P., Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Trans- 
actions on Pattern Analysis and Machine Intelligence 12 (7), 629-639. 

[26] Perona, P., Shiota, T., Malik, J., 1994. Anisotropic Diffusion. Kluwer Academic Press, pp. 73—92. 

[27] Petrila, T., Trif, D., 2005. Basics of Fluid Mechanics and Introduction to Computational Fluid 
Dynamics. Springer. 

[28] Pottmann, H., Wallner, J., Huang, Q., Yang, Y.-L., 2009. Integral invariants for robust geometry 

processing. Computer Aided Geometric Design 26, 37-60. 
[29] Pratt, W. K., 2001. Digital Image Processing, 3rd Edition. John Wiley & Sons. 
[30] Scharr, H., Korkel, S., Jahne, B., 1997. Numerische isotropieoptimierung von fir-filtern mittels 

querglattung. In: Proc. DAGM'97. pp. 367-374. 
[31] Schneider, R., 1993. Convex Bodies: The Brunn-Minkowski Theory. Cambridge University Press. 
[32] Seymour, P. D., Zaslavsky, T., 1984. Averaging sets: A generalization of mean values and spherical 

designs. Advances in Math. 52, 213224. 
[33] Stroud, A. H., 1971. Approximate Calculation of Multiple Integrals. Prentice-Hall. 



20 



[34] Taubin, G., 1995. Estimating the tensor of curvature of a surface from a polyhedral approximation. 

In: Proc. ICCV'95. pp. 902-907. 
[35] Valenti, R., Sebe, N., Gevers, T., 2009. Image saliency by isoccntric curvcdncss and color. In: IEEE 

International Conference on Computer Vision. 
[36] Wardetzky, M., Mathur, S., Kiilbercr, F., E., C, July 2007. Discrete Laplace operators: No free 

lunch. In: Proceedings of Fifth Eurographics Symposium on Geometry Processing. Barcelona, 

Spain, pp. 33-37. 

[37] Weickert, J., 1994. Anisotropic diffusion filters for image processing based quality control. In: Proc. 

Seventh European Conf. on Mathematics in Industry. Teubner-Verlag, pp. 355-362. 
[38] Weickert, J., 1998. Anisotropic Diffusion in Image Processing. Teubner-Verlag, Stuttgart. 
[39] Weickert, J., Scharr, H., 2002. A scheme for coherence-enhancing diffusion filtering with optimized 

rotation invariance. Journal of Visual Communication and Image Representation 13, 103—118. 
[40] Zhong, L., Su, Y., Yeo, S.-Y., Tan, R.-S., Ghista, D. N., Kassab, G., 2009. Left ventricular regional 

wall curvedness and wall stress in patients with ischemic dilated cardiomyopathy. American Journal 

of Physiology - Heart and Circulatory Physiology 296, H573-H584. 



21 



